function [h, n] = inversez(G, Zeros, Poles)
%Author: Paul Ozog
%
%Usage:
% G: Gain.  Must be real.
% Zeros: Vector of zeros.  Complex zeros must be in pairs of complex conjugates.
% Poles: Vector of poles.  Complex poles must be in pairs of complex conjugates.

%Check inputs

if ~(imag(G) == 0)
  disp('Error: G must be real')
  return
end

imagZeros = Zeros(find(imag(Zeros) ~= 0));
imagPoles = Poles(find(imag(Poles) ~= 0));

for i=1:length(imagZeros)
  if ~length((find(imagZeros == conj(imagZeros(i)))))
    disp('Error: Zeros are not in conjugate pairs');
    return;
  end
end

for i=1:length(imagPoles)
  if ~length((find(imagPoles == conj(imagPoles(i)))))
    disp('Error: Poles are not in conjugate pairs');
    return;
  end
end

leftPoles = Poles(abs(Poles) > 1).^-1;
rightPoles = Poles(abs(Poles) <= 1);
leftZeros = Zeros(abs(Zeros) > 1).^-1;
rightZeros = Zeros(abs(Zeros) <= 1);

leftNum   = poly(leftZeros);
leftDenom = poly(leftPoles);
rightNum = poly(rightZeros);
rightDenom = poly(rightPoles);

done = 0;
L = 10;

while ~done

  hL = impz(leftNum, leftDenom,  L+1);
  if ~length(leftPoles)
    hL = fliplr(hL');
  else
    hL = -1*fliplr(hL');
  end

  hR = impz(rightNum,rightDenom, L+1);
  hR = hR';

  hLinv = impz(leftDenom, leftNum, L+1);
  if ~length(leftPoles)
    hLinv = fliplr(hLinv');
  else
    hLinv = -1*fliplr(hLinv');
  end

  hRinv = impz(rightDenom, rightNum, L+1);
  hRinv = hRinv';

  h = G * conv(hL, hR);
  hinv = G^-1 * conv(hLinv, hRinv);
  n = [-L:L];
  
  d = conv(h,hinv);

  d_0 = d(2*L+1)^2;
  Gamma = abs( d_0 / ((sum(d.^2) - d_0)));
  disp(['Curent value of d[n] / (sum(d^2[n] - d^2[0])) : ', num2str(Gamma)])

  if (Gamma > 10^4)
    done = 1;
  end
  
  %if Gamma approaches 0, then there's a problem so just exit the loop
  if (Gamma <= 1)
    done = 1;
  end

  %double L if not done and try again
  if ~done
    L = L*2;
  end
  
end %end while

subplot(1,3,1);
stem(n, h, '.');
title('h[n]');
xlabel('n');
subplot(1,3,2);
stem(n, hinv, '.');
title('h_{inv}[n]');
xlabel('n');
subplot(1,3,3);
stem([-2*L:2*L], d, '.');
title('\delta[n] Approximation');
xlabel('n');
